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Abstract 

Sex chromosome divergence has been documented across phylogenetically diverse species, with amphibians typically having cyto- 
logically nondiverged (" homomorphic") sex chromosomes. With an aim of further characterizing sex chromosome divergence of an 
amphibian, we used "RAD-tags" and Sanger sequencing to examine sex specificity and heterozygosity in the Western clawed frog 
Silurana tropicalis (also known as Xenopus tropicalis). Our findings based on approximately 20 million genotype calls and approxi- 
mately 200 polymerase chain reaction-amplified regions across multiple male and female genomes failed to identify a substantially 
sized genomic region with genotypic hallmarks of sex chromosome divergence, including in regions known to be tightly linked to the 
sex-determining region. We also found that expression and molecular evolution of genes linked to the sex-determining region did not 
differ substantially from genes in other parts of the genome. This suggests that the pseudoautosomal region, where recombination 
occurs, comprises a large portion of the sex chromosomes of 5. tropicalis. These results may in part explain why African clawed frogs 
have such a high incidence of polyploidization, shed light on why amphibians have a high rate of sex chromosome turnover, and raise 
questions about why homomorphic sex chromosomes are so prevalent in amphibians. 

Key words: sex chromosome, pseudoautosomal region, recombination, sex determination, African clawed frogs, Xenopus 
tropicalis. 



Introduction 

Sex can be advantageous because it decouples beneficial from 
deleterious mutations via recombination, which increases the 
variance in fitness effects of linked mutations, and thus the 
efficiency with which natural selection operates. In species 
with genetic sex determination, developmental differences 
between the sexes are initiated by genetic differences be- 
tween the sex chromosomes. In some lineages, the genes 
responsible for triggering sex determination vary, and the 
sex chromosomes (which carry the sex-determining region) 
are routinely reassigned from one to another ancestral pair 
of autosomal chromosomes (Fridolfsson et al. 1 998; Ross et al. 
2009; Evans et al. 2012; Pease and Hahn 2012). Ironically, 



suppression of recombination within a sex-specific region is 
often favored by natural selection, lest a sex-specific, sex-de- 
termining allele loses its sex specificity. 

The origin of sex chromosomes could be initiated by sexual 
antagonism (van Doom and Kirkpatrick 2007), and in many 
species, this is associated with cessation of recombination be- 
tween a portion of the sex chromosomes that makes possible 
unisexual inheritance of a key genomic region that triggers sex 
determination. Cessation of recombination between the sex 
chromosomes can be achieved by reducing or eliminating ho- 
mology (Charlesworth 1991), for example, through point mu- 
tations, inversion, deletion, or insertion of DNA. Strikingly, the 
extent of the nonrecombining region may increase overtime, 
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although this is not necessarily the case (Charlesworth et al. 
2005; Bergero and Charlesworth 2009). The expansion of 
nonrecombining regions may be influenced by the nature of 
evolution in nonrecombining genomic regions, which is 
influenced by Muller's ratchet, background selection, Hill- 
Robertson effects, and genetic hitchhiking of deleterious 
alleles with beneficial mutations (Charlesworth B and 
Charlesworth D 2000). Suppressed recombination between 
sex chromosomes thus has important implications for 
genome evolution, speciation, and adaptation. 

Sex chromosome "degeneration" can be associated with 
sex chromosome divergence resulting from suppressed re- 
combination and involves the loss of coding regions, the ac- 
cumulation of repetitive regions, and structural changes such 
as insertions, deletions, and inversions on the sex-specific 
chromosome (the Y or W). Thus, degenerate sex chromo- 
somes have differences that extend beyond the fundamental 
difference in the presence or absence of a sex-determining 
allele. However, evolutionarily young sex chromosomes start 
out being similar to each other because they originated from 
an essentially identical pair of autosomal chromosomes. In the 
medaka fish, for example, the male-specific region of the Y 
chromosome contains a newly evolved sex-determining locus 
(dmrtlbY) and is only approximately 258,000-bp long (Kondo 
et al. 2006). Similarly, the sex chromosomes of the tiger puf- 
ferfish appear to be distinguished by only one nonsynon- 
ymous substitution (Kamiya et al. 2012). In theory, natural 
selection may drive the expansion of the sex-specific region 
of suppressed recombination (Charlesworth et al. 2005), and 
old sex chromosomes may evolve distinct suites of genes with 
unique, nonhomologous, and/or extensively diverged func- 
tions. In humans, for example, almost all the ancestral genes 
persist on the X chromosome but have been lost on the Y 
chromosome (Skaletsky et al. 2003). The limit of divergence is 
achieved if the sex-specific sex chromosome is lost altogether, 
as occurred in the Ryukyu spiny rat (Kuroiwa et al. 2010). In 
contrast, however, the sex chromosomes of ratite birds and 
boid snakes are old, but each pair is morphologically nondi- 
verged (homomorphic), suggesting that the size of the region 
of suppressed recombination on the sex-specific sex chromo- 
some does not necessarily expand over time (Matsubara et al. 
2006; Tsuda et al. 2007). 

Frog Sex Chromosomes 

Most species of amphibians have homomorphic sex chromo- 
somes (reviewed in Schmid et al. 2010). One possible expla- 
nation for the high incidence of homomorphic sex 
chromosomes in amphibians is that their sex chromosomes 
tend to be young because there has been frequent switching 
of the sex-determining locus during evolution. Consistent with 
this explanation is the inference that sex chromosomes have 
changed ("turned over") approximately 32 times or more 
during evolution based on variation in male (XY) versus 



female (ZW) heterogamy (Evans et al. 2012; Schmid et al. 
2010). Another explanation for widespread homomorphic 
sex chromosomes is that there is periodic recombination 
between the sex chromosomes over most of their length with- 
out changes in the sex-determining locus. In one group of 
hylid frogs, for instance, genomic regions that are tightly 
linked to the sex-determining locus are not substantially 
diverged between males and females, indicating that 
the sex chromosomes of these frogs recombine, at least oc- 
casionally, over most of their length (Stock et al. 201 1). This 
suggests that sex chromosomes of these frogs have large 
"pseudoautosomal" regions where inheritance of genetic in- 
formation resembles autosomal genes and where recombina- 
tion prevents divergence between the sex chromosomes. 
Thus, frequent turnover and recombination clearly both play 
a role in homomorphy of amphibian sex chromosomes, but 
which phenomenon plays the dominant role remains an open 
question. 

The only known amphibian sex-determining gene is called 
DM-W and was discovered in the African clawed frog 
Xenopus laevis (Yoshimoto et al. 2008). DM-W is female 
specific and originated after divergence from the sister 
genus Silurana but before diversification of most or all 
extant species of Xenopus (Bewick et al. 2011). Because 
Silurana tropicalis (also known as X. tropicalis) lacks DM-W, 
sex determination in this species must be triggered by another 
as yet unidentified genetic locus. A high-quality draft genome 
sequence is available for 5. tropicalis that was generated from 
a female (Hellsten et al. 2010), but the sex-specific region of 
this genome has not been characterized. Using amplified frag- 
ment length polymorphisms (AFLPs), Olmstead et al. (2010) 
identified 22 AFLPs linked to the sex-determining locus in the 
"golden" strain of 5. tropicalis (table 1) and proposed that 
females are the heterogametic sex in this strain. Four of 
these 22 AFLPs placed to the distal tip of chromosome/linkage 
group 7 in a linkage map developed by Wells et al. (201 1), also 
represented by scaffold 7 in version 7.1 of the 5. tropicalis 
genome sequence. However, the linkage map contains a 
large (15cM) gap between the most distal two markers 
where the sex-determining region is likely to reside (Wells 
et al. 2011). Many of the other sex-linked AFLPs identified 
by Olmstead et al. (2010) map to other major or small chro- 
mosome/linkage groups (table 1); this is presumably because 
some scaffolds are chimerical (e.g., a portion of scaffold 2 in 
version 7.1 is probably actually derived from S. tropicalis chro- 
mosome 7) and because linkage relationships between some 
small contigs ("orphan scaffolds") and the larger scaffolds 
have not yet been established. Furthermore, additional exper- 
iments with other strains suggest that sex determination may 
occur through the action of multiple alleles at one locus or 
multiple tightly linked genes (Olmstead A, personal 
communication). 

The goal of this study is to further characterize the sex 
chromosomes of 5. tropicalis in terms of the size and level 
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Table 1 

Genomic Regions of Silurana tropicalis That Are Putatively Sex Linked Based on Linkage Study of Olmstead et al. (2010) and Sequencing of 
Chromosome Arm 7p by Seifertova et al. (submitted) 



AFLP 


Recombination 


v4 


v7.1 


7p? 


Portion Sex Linked 


E33.M72.143 


0 


605:241571-241691 


7:4966286-4966166 


Yes 


4435335-5175370 


E33.M81.275 


0 


494:27646-27898 


No hits 


— 


NA 


E33.M90.327 


0 


No hits 


211:76840-76535 


Yes 


ALL 


E38.M93.218 


0 


953:138210-138402 


278:109297-109490 


No 


ALL 


No name 


0 


494:31633-32115 


No hits 


— 


NA 


No name 


0 


494:27541-27902 


No hits 


— 


NA 


No name 


0 


379:817889-818000 


2:149826489-149826600 


Yes 


149787496-150105127 


No name 


0 


736:292586-293067 


78:248343-247864 


Yes 


ALL 


No name 


0.3 


605:245039-245621 


7:4963243^1962661 


Yes 


4435335-5175370 


No name 


0 


605:116800-117215 


22:1040592-1040177 


Yes 


ALL 


E40.M52.572 


0.4 


859:57522-58049 


7:3195527-3196069 


Yes 


76939-3370464 


E33.M61.177 


0.5 


1778:6156-6312; 1778:9971-9815 


144:136541-136385 


Yes 


ALL 


E33.M61.177 


0.5 


1778:6156-6312; 1778:9971-9815 


662:22621-22465 


No 


NONE because 114 is on 7p 


No name 


0.5 


810:261995-262744 


94:264236-264985 


Yes 


ALL 


No name 


0.5 


810:261995-262744 


94:263563-263460 


Yes 


ALL 


E32.M94.406 


1.9 


Multiple hits 


22:58856-59241 


Yes 


ALL 


E37.M52.423 


1.9 


810:325559-325959 


144:95531-95931 


Yes 


ALL 


E37.M52.423 


1.9 


1151:130316-130719 


144:95931-95531 


Yes 


ALL 


E41.M83.506 


1.9 


810:276803-277288 


94:279867-279382 


Yes 


ALL 


E32.M35.552 


2.6 


Multiple hits 


No hits 




NA 


E37.M60.232 


2.6 


6092: 2392-2601 


7931:611-828 


Yes 


ALL 


E32.M59.335 


2.9 


735:292141-292443 


7:7903155-7902853 


Yes 


6687308-9940823 



Note. — NA, not applicable. AFLP refers to the name of the AFLP from Olmstead et al. (2010) if provided. Recombination refers to the recombination rate with the sex- 
determining locus from that study. Scaffold and position of AFLPs are provided for genome assembly version 4.0 (v4) and 7.1 (v7.1). Sex-linked portions that were included in 
categories in tables 2 and 3 based on the level of recombination (Portion sex linked, with NA meaning not applicable) either refer to base pair positions of a contig within a 
larger scaffold that is not interrupted by unknown sequence or the entire scaffold was assumed to be sex linked (ALL). 



of divergence of the sex-specific region, and to compare mo- 
lecular evolution and expression of sex-linked and nonsex- 
linked genes. To this end, we used restriction-site associated 
DNA (RAD) tags (Baird et al. 2008), a reduced representation 
next-generation sequencing approach, to genotype millions of 
homologous nucleotide positions in male and female individ- 
uals including positions that are monomorphic in both sexes, 
polymorphic in one or both sexes, and positions in which a 
genotype inference (i.e., homozygous or heterozygous) was 
only possible in one sex due either to sex specificity of the 
genotyped position or differences in coverage of that position 
between the sexes. The RAD tag approach produces se- 
quences of thousands of small regions that are adjacent to a 
rare cutting restriction enzyme site. Because the sequenced 
portions of the genome are associated with restriction enzyme 
sites, many homologous sequences are obtained from multi- 
ple individuals. Missing data among individuals can arise in 
unusual cases where mutation generates polymorphism in 
the presence or absence of the restriction enzyme sites or 
because of variation among individuals in the depth of se- 
quencing coverage for a particular region. Our analysis incor- 
porated information on sex-linked regions from Olmstead 
et al. (2010), information from a laser-dissected chromosome 
arm 7p from a male individual (Seifertova et al., submitted), 



which is linked to the sex-linked region identified by Olmstead 
et al. (2010), and the most recent genome assembly (version 
7.1, reference accession PRJNA12348). This study thus pro- 
vides, for the first time, a comprehensive perspective on the 
extent of sex chromosome divergence in this species by eval- 
uating the distribution of homozygous and heterozygous ge- 
notypes, molecular evolution, and gene expression of sex 
chromosomes in the context of the rest of the genome. 

Materials and Methods 

Four female and four male S. tropicalis individuals were ob- 
tained from Xenopus Express (Brooksville, FL). Sex was con- 
firmed by dissection and species assignment achieved by 
comparing between 809 and 812 bp of mitochondrial DNA 
sequence from a portion of the 1 6S gene from each sample to 
homologous sequence data from all other known species of 
African clawed frog (Evans et al. 201 1). We performed a phy- 
logenetic analysis on these 8 sequences, 27 sequences from 
individuals used in the polymerase chain reaction (PCR) screen 
detailed earlier, all Silurana sequences from Evans etal. (2004), 
6 5. tropicalis samples from Ghana (obtained from tissue 
archive at the Burke Museum, University of Washington, 
accession numbers UWBM5957-8, UWBM5961-63, and 
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UWBM5969), and sequences from 6 individuals from the 
"golden" strain used by Olmstead et al. (2010) that were 
provided by Richard Harland. We used an X. laevis sequence 
from South Africa as an outgroup in this analysis, and the total 
alignment length was 81 7 bp. Model selection for phyloge- 
netic analysis was accomplished using MrModeltest2 
(Nylander 2004). Phylogenetic analysis was performed with 
MrBayes version 3.1.2 (Huelsenbeck and Ronquist 2001) 
using the best-fit model based on the Akaike Information 
Criterion, with two independent Markov chain Monte Carlo 
(MCMC) runs, each for 2,000,000 generations. Convergence 
of the MCMC runs on the posterior distribution was assessed 
by inspecting parameter trends and effective sample sizes 
using Tracer version 1.5 (Rambaut and Drummond 2007). 
Based on these analyses, a burn-in of 500,000 generations 
was discarded before constructing a consensus tree with 
MrBayes. 

Genomic DNA was extracted from liver using QIAGEN 
DNeasy kit, purified using QIAGEN's spin purification protocol, 
and RAD tag library preparation performed by Floragenex, Inc 
(Eugene, OR). For each individual, two libraries were gener- 
ated — one used the restriction enzyme Sbfl and another used 
Notl. The RAD tag libraries were multiplexed on three lllumina 
flow cells using individual barcodes, and lllumina sequencing 
was performed at the University of Oregon. These data have 
been deposited in GenBank (accession number SRP022004). 

lllumina sequence reads were sorted by barcode with 
RADtools v1 .2.4 using the "fuzzy_MID" option, which assigns 
reads with barcode errors to the nearest barcode (Baxter et al. 
201 1). Data from each individual were independently aligned 
to the S. tropicalis version 7.1 genome using bwa-0.6.2 (Li and 
Durbin 2009) and samtools.0.1 .18 (Li et al. 2009). The 
"MarkDuplicates" function in picard (http://picard.source 
forge.net) was used to mark putative PCR-amplified dupli- 
cates, which were then excluded from the genotyping analysis 
with an aim of minimizing genotyping error. The Genome 
Analysis Toolkit (GATK) version 2.2-15 was then used to re- 
align indels using the "RealignerTargetCreator" and 
"IndelRealigner" functions (McKenna et al. 2010; DePristo 
et al. 2011). The "FixMatelnformation" function of picard 
was then used to adjust mate pair alignments. 

Following "Best Practices" guidelines on the GATK website 
and forum (http://gatkforums.broadinstitute.org/) for analysis 
of genomes that lack known single-nucleotide polymorphisms 
(SNPs), the "UnifiedGenotyper," "BaseRecalibrator," and 
"PrintReads" functions of GATK were used to iteratively ge- 
notype, recalibrate base quality scores, and generate new 
input (bam) files, using the genotype files generated from 
"UnifiedGenotyper" as known polymorphic positions to be 
ignored for base recalibration in each iteration. Convergence 
was reached by the fifth iteration, in that variable positions 
recovered from this analysis were 99.8% identical to those 
from the fourth iteration. The "VariantFiltration" and 
"SelectVariants" functions of GATK were then used to 



identify and exclude genotyped positions that 1) were 
within 1 0 bp of an insertion/deletion, 2) had a Phred genotype 
quality score (Ewing and Green 1998) of less than 30, which 
means that we removed positions that had a probability of 
error of greater than 0.001 , or 3) had more than one-tenth of 
the reads mapping equally well to another position and where 
there were at least four of these reads. 

The S. tropicalis genome assembly 7.1 consists of 7,730 
scaffolds aggregated from 55,234 contigs connected by 
"N"s within each scaffold. The total number of bases is 
1,437,594,934, of which 5% (n = 71, 599,926) are "N"s. 
This assembly includes 14 large "super scaffolds" that were 
assembled using meiotic map, synteny, and cytological data, 
corresponding to the 10 haploid chromosomes, with some 
chromosomes being represented by multiple scaffolds (3a 
and 3b; 5a and 5b; and 8a, 8b, and 8c). The rest of the scaf- 
folds are "orphan scaffolds" whose chromosomal locations 
are not yet known. We divided the genomic regions into five 
mutually exclusive groups based on 1) the inferred level of 
recombination with the sex-determining region by Olmstead 
et al. (2010), 2) the linkage groups in the genome assembly 
7.1 (table 1), and 3) the results of the lllumina sequencing of 
the dissected petite arm of chromosome 7 (Seifertova et al., 
submitted). The first of the five groups ("completely sex 
linked") included contigs from assembly 7.1 that contain re- 
gions that had no recombination (0%) with the sex-determin- 
ing region in Olmstead et al. (2010). This means that 
recombination between an AFLP polymorphism and the sex- 
determining region was not observed in any of 300 individuals 
assayed by Olmstead et al. (2010). The second group ("par- 
tially sex-linked") included contigs from assembly 7.1 that 
contain regions that had a recombination rate of more than 
0% and less than 3.0% in Olmstead et al. (2010). The third 
group ("chromosome 7p") contained sections of scaffolds in 
assembly 7.1 that are located on chromosome 7p according 
to Seifertova et al. (submitted), and not in the "completely sex 
linked" or "partially linked" categories. The fourth group 
("non-7p chromosomes") contained the remaining sections 
on the chromosome-scale scaffolds in Assembly 7.1, including 
the portion of scaffold 7 that did not map to chromosome 7p. 
The fifth group ("other orphans") contained orphan scaffolds 
in assembly 7.1 that 1) have not been linked to a chromo- 
some, 2) have no evidence of sex linkage according to 
Olmstead et al. (2010), and 3) did not map to chromosome 
arm 7p according to Seifertova et al. (submitted). More spe- 
cific information on the scaffold or scaffold portions in each of 
these groups is provided in table 1 . 

Genome-Wide Distribution of Genotypes in Female and 
Male S. tropicalis 

For genomic regions in each of the six categories described 
earlier, we tabulated genotype patterns for three scenarios 
(fig. 1) in 500,000-bp windows across the S. tropicalis 
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Scenario la: Scenario 2a: Scenario 3a: 

W/Z divergence Z-specificity W-specificity 

z Z 

w 

9 z= z 

Expectation for O >> Cf O =0 Q-Oand 

heterozygosity: * V no male tags 

Scenario 1b: Scenario 2b: Scenario 3b: 

X/Y divergence X-specificity Y-specificity 

o x x 

V X x 

rf y y 

U x x 

Expectation for O << Cf Cf=° C? = 0and 

heterozygosity: V u no female tags 



Fig. 1. — Genotypic scenarios for sex-linked regions. Expectations for 
heterozygosity depend on which sex is heterogametic and the region of 
the sex chromosome (pseudoautosomal region vs. sex-determining 
region). Female heterogamy is potentially associated with female-biased 
heterozygosity (Scenario 1a), male-only heterozygosity (Scenario 2a), or 
female-only homozygosity with no male genotypes (Scenario 3c). 
Corresponding scenarios (Scenarios 1b, 2b, and 3b) apply to the opposite 
sex for male heterogamy. 



genome; smaller windows were examined at the ends of scaf- 
folds or when a scaffold was smaller than 500,000 bp. 
Genotype patterns in each sex (i.e., the distribution of homo- 
zygous or heterozygous positions) are relevant to sex chromo- 
some evolution in the following ways. First, divergence 
between the sex chromosomes due to suppressed recombi- 
nation generates positions that are either heterozygous in all 
females and no males (for a ZW sex-determining system) or 
heterozygous in all males and no females (for an XY sex- 
determining system). We call these patterns "Scenario 1a" 
and "Scenario 1b," respectively (fig. 1). We note that in 
"Scenario 1a" regions, some positions can also be heterozy- 
gous in males due to polymorphism on the Z chromosome, 
and in "Scenario 1b," some positions can also be heterozy- 
gous in females due to polymorphism on the X chromosome. 
In any case, in genomic regions consistent with Scenario 1 , 
heterozygosity observed in all samples from one sex is ex- 
pected to exceed heterozygosity observed in all samples 
from the other sex. We, therefore, searched for regions with 
heterozygosity present in all females or in all males. For both of 
these statistics, we ignored positions that are heterozygous in 
all genotyped individuals. To account for variation in coverage 
in males and females, for each window, we divided these 
counts by the total number of positions in each window for 
which genotype calls were made in at least one female and at 
least one male. 



Another genotypic scenario for sex chromosomes is that a 
genomic region may be present only on the Z chromosome 
(with female heterogamy) or only on the X chromosome (with 
male heterogamy) (Scenarios 2a and 2b; fig. 1). No counter- 
part exists on the W chromosome (or Y chromosome) due to 
deletion, insertion, or divergence. To detect such a genomic 
region, we searched for regions with heterozygous positions 
present in one sex but not the other. For such positions, we 
required a genotype call in at least one individual of each sex 
but heterozygous calls to be present in only one sex. To ac- 
count for variation in coverage in males and females, for each 
window we divided these counts by the total number of po- 
sitions in each window for which genotype calls were made in 
at least one female and at least one male. 

A third genotypic scenario for sex chromosomes is that a 
genomic region may be present only on the W chromosome 
or only on the Y chromosome (Scenarios 3a and 3b, fig 1). 
Thus, we searched for positions that had genotype calls only in 
females (or only in males) and that are all homozygous. To 
account for variation in coverage, we standardize the counts 
in each window by the sum of the number of positions in each 
window for which genotype data are available for 1) at least 
one female and at least one male, 2) at least one female but 
no males, and 3) at least one male but no females. Thus, by 
evaluating these three genotype scenarios in genomic win- 
dows across the 5. tropicalis genome assembly, we attempted 
to identify genomic windows that either had significantly 
more heterozygous positions in one sex (Scenario 1), that 
had heterozygous positions only in one sex (Scenario 2), or 
that had homozygous positions in only one sex and no ho- 
mologous genotypes in the other (Scenario 3). Higher values 
for each ratio are suggestive of genotype patterns character- 
istic of diverged sex chromosomes. 

Expression and Molecular Evolution 

As described in Chain et al. (2011), we estimated gene ex- 
pression levels based on sequences across 26 expressed 
sequence tag (EST) libraries from the following tissues or de- 
velopmental stages: egg, gastrula, neurula, embryo, tailbud, 
tadpole, metamorphosis, adipose tissue, bone, brain, head, 
heart, intestine, kidney, limb, liver, lung, ovary, oviduct, skel- 
etal muscle, skin, spleen, stomach, tail, testis, and thymus. We 
summarized patterns of gene expression across EST libraries 
using the nonindependent "total," "intensity," and "even- 
ness" statistics described in Chain et al. (2011). The "total" 
expression of a gene (7) is the proportion of times that a gene 
was sequenced in each EST library (/_,) summed across all li- 
braries (T= J^LI). The "intensity" of expression (I) is the mean 
expression level from the perspective of a gene and is calcu- 
lated following this equation: l=Y^U 2 ^HU "Evenness" of 
expression (£) can be thought of as the "effective number" 
of tissues in which a gene is expressed and is calculated 
following this equation: E = Til. 
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For a subset of the sex-linked and nonsex-linked genes, we 
also calculated the rate ratio of nonsynonymous to synony- 
mous substitutions per site (d/V/dS) along the S. tropicalis lin- 
eage using PAML version 4.5 (Yang 1997). This ratio was 
calculated using a maximum likelihood model that individually 
estimates d/V/d5 for each branch in a phylogeny, following 
Chain et al. (201 1). Our phylogeny was estimated from se- 
quences from 5. tropicalis, X. laevis, and using sequences from 
another pipid frog (Pipa carvalhoi or Hymenochirus curtipes) as 
an outgroup. To avoid undefined values, we added 0.02 to all 
dS values before calculating d/\//d5, following Chain et al. 
(201 1). We made this adjustment a priori by looking only at 
dS values, to make better use of the data. Because extreme 
values for d/V and dS were occasionally estimated, we ex- 
cluded from the analysis genes with an estimated d/V or dS 
value above 2, and any genes whose available data comprised 
less than 1 00 synonymous positions. 

We used a one-sided permutations to test whether the 
expression and molecular evolutionary statistics differed be- 
tween genes that either (a) were or (b) were not on the same 
chromosome as the sex-determining locus. The permutations 
randomly divided the set of (a + b) values into two groups of 
size a and b and then calculated the difference between the 
averages of each group. We repeated this 1,000 times to 
generate a distribution for the null hypothesis that the 
values were drawn from the same underlying distribution, 
and then compared this with the observed differences, 
which is the test statistic of each test. A significant difference 
was inferred if the observed difference was greater than 95% 
of the differences from the permutations. Because these tests 
are one sided, the operands of the test statistic (i.e., the min- 
uend and subtrahend of each difference) were defined ac- 
cording to specific expectations for sex chromosome 
degeneration discussed later. 

Results 

Mitochondrial DNA Variation within S. tropicalis, 
Including the "Golden" Strain 

We analyzed phylogenetic relationships among approximately 
810-bp region of mitochondrial DNA from the commercially 
obtained S. tropicalis individuals that we used for RAD tags 
and PC R screens, six individuals from the golden strain used by 
Olmstead et al. (2010), and several other wild-caught 
5. tropicalis individuals and individuals from other Silurana spe- 
cies. An identical mitochondrial DNA sequence was obtained 
from the six golden strain individuals, one of the samples we 
used for RAD tag sequencing (a female), 20 of the samples we 
used for PCR screens (9 females and 1 1 males), and one indi- 
vidual sampled from Nigeria. Mitochondrial sequences from 
five samples used in the RAD tag sequencing (2 females and 3 
males) were identical to each other and differed from the 
golden strain sequence by one nucleotide substitution. 



Mitochondrial sequences from two other samples used in 
RAD tag sequencing (1 female and 1 male) and seven samples 
used in PCR screens (4 females and 3 males) were identical to 
each other and differed from the golden strain mitochondrial 
sequence by a different single-nucleotide substitution than the 
previously mentioned sequence present in five of the RAD tag 
samples. Mitochondrial sequences from another sample from 
Nigeria differed from the golden strain mitochondrial se- 
quence by two nucleotide substitutions. Phylogenetic analysis 
of these and other sequences indicates that the commercially 
obtained S. tropicalis samples used in this study form a well- 
supported clade that includes two sequences from Nigeria and 
the six sequences from the golden strain of 5. tropicalis (fig. 2). 
This clade is possibly common in individuals east of the 
Dahomey Gap, a savannah corridor that interrupts the West 
African rain forest (Salzmann and Hoelzmann 2005). 

Reduced Representation Genome-Wide Genotyping from 
RAD Tags 

We used a reduced representation genome sequencing ap- 
proach called "RAD tags" to sequence many small but ho- 
mologous portions of the 5. tropicalis genome in four female 
and four male individuals. An average of 9,696,525 lllumina 
reads were mapped in each individual, with the average 
number of reads mapped per female or per male being 
9,445,507 and 9,947,543 reads, respectively. After excluding 
positions in the reference sequence with no data, within an 
individual the average depth of coverage was 18.4 reads per 
position. Genotypes were called for a total of 19,624,824 
positions, and 193,199 SNPs (0.98%) were detected. For 
each position at which at least one genotype was called, an 
average of 7.14 out of 8 individuals were genotyped. 

If S. tropicalis has a large female-specific genomic region on 
the W chromosome, we expected a higher proportion of the 
lllumina reads from females to map to the genome assembly 
because this assembly was generated from a female individ- 
ual. Contrary to this expectation, a slightly higher proportion 
of reads from males (average per male individual 89.2%, 
range: 87.8-90.8%) than from females (average per female 
individual 87.9%; range: 85.0-90.8%) mapped to this 
genome assembly, arguing against there being a large 
female-specific region in the 5. tropicalis genome. Another 
indication of a large female-specific genomic region on the 
W chromosome would be a substantially higher number of 
positions genotyped in females than in males. Out of a total of 
19,624,824 positions that were genotyped with high confi- 
dence in at least one individual, slightly more genotypes were 
recovered in females than in males: 912,738 (4.7%) positions 
were genotyped only in one or more females, and 462,792 
(2.4%) positions genotyped only in one or more males. 
However, in the 10 largest scaffolds, the number of genotype 
calls in at least one female was consistently 1 .1-3.3% higher 
than the number of genotype calls in at least one male, with 
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Fig. 2. — Phylogenetic analysis of mitochondrial DNA sequences suggest that the "golden" strain used by Olmstead et al. (2010) and samples used in this 
study (8 for RAD tag analysis and 27 others for PCR assays) originate from Nigeria. Nodes with >95% posterior probability are indicated with a black circle. 
Species names, including those undescribed, follow Evans et al. (2004). Abbreviated country names include the Republic of the Congo (R. Congo), the 
Democratic Republic of the Congo (DRC), and Equatorial Guinea (EG). The scale bar refers to the number of substitutions per site, gray areas on the map 
indicate the distribution of tropical forest in West Africa, and the dotted line indicates the approximate distribution of Silurana tropicalis inferred by Tinsley 
etal. (1996). 



scaffold 7 having 2.5% more genotype calls in females than 
males. This suggests that the higher number of unique geno- 
type calls in females is primarily a technical artifact related to 
differences in coverage among individuals in the RAD tag li- 
braries. The RAD tag data did not provide high-quality geno- 
types from any positions on 5,721 scaffolds, which together 
comprise 36,133,437 bp (-2.1 % of the genome). 

Genome-Wide Genotype Patterns and Nucleotide 
Diversity Similar in Males and Females 

We searched 500,000-bp windows for various genotypic pat- 
terns consistent with sex chromosome divergence expected 
under female or male heterogamy (fig. 1). In general, this 
effort failed to identify any regions with a pronounced geno- 
typic signature of sex chromosome divergence expected by 
female heterogamy (table 2). One exception was a significant 
excess of windows with female-only homozygous genotypes 
(Scenario 3a) in orphan scaffolds, but we suspect this was an 
artifact related to the broader coverage in females. Most no- 
tably, portions of linkage groups 2 and 7 that were catego- 
rized as "partially sex-linked" and "completely sex-linked" to 
the sex-determining region based on Olmstead et al. (2010) 



did not exhibit a genotypic pattern consistent with degenerate 
sex chromosomes based on the RAD tag genotypes. 

Considerable caution is needed in the interpretation of the 
average genotype frequencies in genomic windows for the 
"other orphans" category because in many cases the scaffold 
is smaller than the window size (500,000 bp), and the result- 
ing truncated genomic windows are therefore expected to 
have an increased variance in the frequency of various geno- 
typic patterns. Additionally, average genotype frequencies in 
these genomic windows could fail to detect small scaffolds 
that have genotypic patterns consistent with sex chromosome 
divergence. For example, "other orphans" had higher than 
expected values for Scenarios 2b and 3b, which are consistent 
with male heterogamy, and this is probably related to the 
small size of these scaffolds and consequent increase in the 
sex-specific genotypes in truncated windows for these 
scaffolds. 

Additional insights are gained by examining nucleotide di- 
versity in each sex within 500,000-bp windows. If a portion of 
the sex chromosomes is substantially diverged, we expected 
much higher average nucleotide diversity per site in one sex 
(females for female heterogamy) in genomic windows span- 
ning this diverged region. However, average nucleotide 
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diversity per site is essentially identical in males and females 
throughout these scaffolds, including chromosome arm 7p, 
which is linked to the sex-determining region (fig. 3). To ex- 
plore the possibility that there could be variation within the 
RAD tag samples in sex chromosome divergence that corre- 
sponds with the three mitochondrial DNA haplotype groups 
detailed earlier, we explored nucleotide diversity in male and 
female individuals from each group. This analysis also did not 
identify a pronounced signature of sex chromosome diver- 
gence (supplementary figs. S1 and S2, Supplementary 
Material online). 



Genotype Patterns Based on Sanger Sequencing 

We amplified 65 portions (amplicons) of genomic regions 
identified by Olmstead et al. (2010) to be linked to the sex- 
determining region, including 18 and 46 amplicons from 
"completely sex-linked" and "partially sex-linked" regions, 
respectively (table 1, supplementary table S1 and fig. S3, 
Supplementary Material online). None had female-specific 
amplifications, allowing us to dismiss Scenario 3 (fig. 1) for 
all these regions. We sequenced 45 of these amplifications in 
multiple male and female individuals. SNPs or insertion/dele- 
tion polymorphisms were shared between males and females 
in at least one amplicon for essentially all scaffolds (no poly- 
morphism was observed in amplicons from scaffold 144). This 
suggests that Scenario 1 is unlikely for these regions, with the 
caveat being that a heterozygous position could arise in both 
sexes in a region consistent with Scenario 1 through conver- 
gent evolution on the W and Z. 

We also used PCR to examine an additional 173 regions 
that exhibited signs of sex linkage based on our analyses of the 
RAD tag data, including regions of chromosome 7p and else- 
where as detailed in supplementary figure S3 and table S1, 
Supplementary Material online. None had sex-specific ampli- 
fications, allowing us to dismiss Scenario 3 for all these re- 
gions. We sequenced 94 of these amplifications from male 
and female individuals. Thirty of these were not polymorphic 
in any of the individuals we sequenced. Forty-eight had poly- 
morphisms shared between males and females, allowing us to 
conclude that Scenario 1 is unlikely for these regions. Five had 
polymorphisms in both sexes with none being shared across 
sexes. Ten had polymorphisms only in females and two had 
polymorphisms only in males. 

Two amplifications were of particular interest. An amplifi- 
cation on scaffold 7 that spanned positions 10,128,301- 
10,129,920 was highly polymorphic in females but not 
males, although no polymorphism was fixed in females (out 
of seven females and three males sequenced; supplementary 
table S1, Supplementary Material online). This region failed to 
amplify in four females and three males. Another amplifica- 
tion, which targeted a region on scaffold 163 had 31 poly- 
morphisms in three females but only one polymorphism in 
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Fig. 3. — Average nucleotide diversity per site (n) in 500,000-bp windows is similar in males (blue) and females (red) throughout much of the Silurana 
tropicalis genome. Plots are labeled with numbers referring to scaffolds 1-1 0, which collectively comprise approximately 75% of the genome. A gray bar on 
scaffold/chromosome 7 indicates the petite arm based on the linkage map of Wells et al. (201 1 ), which carries the sex-determining region in the 5. tropicalis 
golden strain (Olmstead et al. 2010). 



three males, but essentially all the female polymorphisms were 
present in only one individual. 

Gene Expression and Molecular Evolution 

Expression was detected in a total of 37,790 transcripts in at 
least one of the 26 EST libraries we surveyed (table 3). On the 
basis of studies of recently diverged neosex chromosomes in 
fruit flies (Drosophila) (reviewed in Bachtrog 2013), we ex- 
pected expression of genes situated near the sex-determining 
locus to be expressed 1) at a lower total level, 2) higher inten- 
sity, 3) lower evenness, and to have 4) higher d/V/dS compared 
with genes in other parts of the genome (i.e., the "non 7p 
chromosomes"). For the most part, these expectations were 
not met for genes that were demonstrably very close to the 
sex-determining region, with the one exception that the in- 
tensity of "completely sex-linked" genes were individually sig- 
nificantly higher than the "non 7p chromosomes" (P< 0.05, 
table 3). Evenness of "other orphans" was also significantly 
lower than "non 7p chromosomes" as was total expression of 
"other orphans." d/V/dS was significantly higher only in 



"chromosome 7p" compared with "non-7p chromosomes" 
but the magnitude of this difference was small. No expression 
data were recovered from genes on Scaffold 22, which is 
tightly linked to the sex-determining region (Olmstead et al. 
2010), even though it was 1,1 56,260-bp long. We examined 
this scaffold using Xenbase (Bowes et al. 2009) and found that 
it contained a cluster of olfactory receptors, which (not sur- 
prisingly) were not highly expressed in any of the EST libraries 
we examined. 

Discussion 

To explore sex chromosome divergence in an amphibian, we 
used genotype calls from approximately 20 million positions, 
information about sex linkage, EST databases, and molecular 
evolutionary analyses to further characterize the sex chromo- 
somes of the Western tropical frog S. tropicalis. Phylogenetic 
analysis of mitochondrial DNA sequences suggests our sam- 
ples originated in Nigeria, which is also the source of the 
female individual from which the genome sequence was gen- 
erated (Hellsten et al. 2010). Additionally, our analysis also 
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Table 3 



Average Expression and Molecular Evolution Statistics for Silurana tropicalis Genes in Five Genomic Categories 



Region 


Number of Genes 
(Expression) 


Total 


Intensity 


Evenness 


Number of 
Genes (d/V/dS) 


d/v/dS 


"Non-7p chromosomes" 


35,135 


0.00065 


0.00015 


3.19655 


9183 


0.2702 


"Chromosome 7p" 


2,114 


0.00079 


0.00019 


3.17054 


546 


0.2850* 


"Other orphans" 


260 


0.00040* 


0.00012 


2.59443* 


55 


0.2670 


"Partially sex-linked" 


246 


0.00090 


0.00018 


3.37658 


71 


0.2751 


"Completely sex-linked" 


35 


0.00101 


0.00049* 


2.56020 


8 


0.2615 



Note. — See Materials and Methods for description of statistics. 

*Values that are individually significantly different from the "Non-7p chromosomes" (P<0.05, one-sided permutation tests). 



suggests that the golden strain analyzed by Olmstead et al. 
(2010) is from Nigeria. 

Known sequences in the S. tropicalis genome sequence 
assembly version 7.1 comprise approximately 80.4% of the 
approximately 1 .7 Gbp genome, and scaffolds, including 
"N"s, comprise approximately 84.5% of the genome. Thus, 
the RAD tag data could not be compared with 1 5-20% of the 
genome because of gaps in the genome sequence. Because of 
variation in coverage, high confidence genotype calls were not 
made on scaffolds that together comprise an additional 2.1 % 
of the genome. Thus, in this study, we lack information from a 
nontrivial portion of this genome. 

Mindful of these substantial gaps in genome sequence and 
the uncertainty in linkage relationships among many unas- 
sembled (orphan) scaffolds, we leveraged information from 
a targeted sequencing effort of chromosome arm 7p and also 
the linkage analysis by Olmstead et al. (2010) to guide our 
analysis. The dearth of genotypic patterns consistent with di- 
vergent sex chromosomes, and particularly patterns that are 
consistent with female heterogamy (table 2), and the similar 
level of pairwise nucleotide diversity in males and females 
throughout the petite arm of chromosome 7 (fig. 3) argues 
strongly against there being a large sex-specific region of the 
S. tropicalis chromosomes. This inference is consistent with the 
findings of Uno et al. (2008) who detected no sex differences 
in C-banded heterochromatin in 5. tropicalis. 

On the basis of studies of fruit flies (reviewed in Bachtrog 
201 3), we expected genes linked to the sex-determining locus 
to potentially exhibit lower total expression and higher speci- 
ficity (i.e., higher intensity and lower evenness as defined in 
Materials and Methods section). We also expected molecular 
evolution of these genes to be consistent with relaxed purify- 
ing selection. However, on the basis of a small sample size, we 
only observed a significant increased expression intensity of 
"completely linked" genes compared with the rest of the 
genome, with none of these expectations met in "partially 
sex-linked" genes (table 3). Some of these expectations 
were also met in orphan scaffolds, which have undetermined 
linkage relationships with respect to the sex-determining 
locus, and regions of chromosome arm 7p. It is not clear 



that these latter observations are related in any way to linkage 
to the sex-determining region. 

Caveats exist in our interpretation of these data. First, non- 
recombining portions of the genome tend to accumulate re- 
petitive sequences that can be difficult to sequence and map. 
For this reason, the sex-specific portion of the S. tropicalis 
genome may be under-represented in the current genome 
assembly and/or our mapped lllumina reads. Second, it 
is conceivable that there is polymorphism in the sex- 
determining mechanism (Olmstead A, personal communica- 
tion). Polymorphism in genetic sex determination could occur 
at a single locus wherein multiple, differently functioned sex- 
determining alleles are segregating at a single locus, which 
have distinct and not necessarily transitive dominance relation- 
ships. Polymorphism in genetic sex determination could also 
occur at multiple loci distributed on the same or different 
chromosomes. Sex determination in zebrafish, for example, 
appears to be orchestrated by genes on different chromo- 
somes (Anderson et al. 2012). Genotypic patterns expected 
with these types of polymorphisms are unclear and could in- 
clude a dearth or absence of pronounced sex chromosome 
divergence. A third caveat to our conclusions is that polymor- 
phism among females could also potentially exist in the extent 
of divergence between the W and Z chromosomes. Under this 
scenario, it is conceivable that there could be variation among 
populations in the extent of recombination along the sex chro- 
mosomes and consequently the extent and magnitude of di- 
vergence between the sex chromosomes. Further exploration 
of these possibilities will be assisted by the identification of the 
sex-determining locus in S. tropicalis, the completion of high- 
quality sequencing and assembly of sex-linked regions, and 
the exploration of variation within and among populations 
in sex determination and sex chromosome evolution. 

Polyploidization, Dosage Compensation, and Sex 
Chromosome Turnover 

Within a species, the propensity to undergo genome duplica- 
tion and sex chromosome evolution is potentially interrelated. 
For example, polyploidization might be less common in spe- 
cies with divergent sex chromosomes where one has 
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degenerated because, after duplication, the degenerate an- 
cestral sex chromosome would segregate as a new autosomal 
chromosome, and the resulting homozygous null genotypes 
could be detrimental (Evans et al. 2012). Sex chromosome 
degeneration also creates imbalances in allelic copy number 
between the sexes, which can lead to the evolution of dosage 
compensation — a factor that is also potentially relevant to 
polyploid speciation (Orr 1990). Dosage compensation is a 
process that equalizes expression levels in each sex of a 
gene that has a different number of alleles in each sex. This 
could evolve in a species with female heterogamy, for exam- 
ple, through inactivation of one of the Z alleles in males or 
through upregulation of the Z allele in females. Orr (1990) 
proposed that dosage compensation in species with a degen- 
erate sex chromosome could act as a barrier to polyploid spe- 
ciation because dosage compensation would be disrupted 
when a newly formed triploid individual backcrosses with a 
diploid parental individual. Our analyses suggest that the sex- 
specific region of S. tropicalis is small, that sex chromosome 
divergence is minimal, and therefore that dosage compensa- 
tion associated with degeneration of the sex-specific sex chro- 
mosome would have evolved in very few genes or none at all. 
Together these features of the sex chromosomes may have 
facilitated (or at least not impeded) polyploidization in 
Silurana, which occurred at least once (reviewed in Evans 
2008). Interestingly, the sister genus Xenopus has a newly 
evolved sex-determining gene called DM-W (Yoshimoto 
et al. 2008; Bewick et al. 2011). Species in this group also 
probably have minimally diverged sex chromosomes and 
have undergone polyploid speciation multiple times (Evans 
2008). Clearly, however, this is not the only consideration in 
the ability of species to tolerate polyploidization because many 
amphibian groups that have homomorphic sex chromosomes 
lack polyploid species. 

The extent of sex chromosome degeneration is also rele- 
vant to the chances a species experiences future sex chromo- 
some turnover — a change in which pair of chromosomes 
carries the trigger for sex determination (Charlesworth and 
Mank 2010). If sex chromosome turnover occurs in a species 
with a diverged and degenerate sex chromosome, the ances- 
tral degenerate chromosome could segregate autosomally, 
and some individuals could inherit two copies and be homo- 
zygous for degenerate alleles (Charlesworth and Mank 2010). 
Thus, sex chromosome turnover may be more likely in species 
that have sex chromosomes that are not substantially 
degenerated. 

If sex chromosome turnover were common, this could 
maintain homomorphy of sex chromosomes. Recent work 
on sex chromosomes in African clawed frogs has established 
nonhomology between the sex chromosomes of X. laevis and 
5. tropicalis (Uno et al., submitted) and the recent appearance 
of a novel sex-determining locus in X. laevis (Yoshimoto et al. 
2008; Bewick et al. 201 1 ). Thus, it appears that the origin of a 
new sex-determining gene in X. laevis was associated with a 



reassignment of sex chromosomes without necessarily involv- 
ing a change in heterogamy, and that recent sex chromosome 
turnover can account for sex chromosome homomorphy in 
this species (Tymowska 1991). In other species, including 
S. tropicalis, it is also possible that the sex-determining mech- 
anism of nondiverged sex chromosomes could be old, but that 
divergence is prevented by periodic recombination, which 
possibly could be facilitated by breeding individuals that are 
phenotypically sex reversed (the "fountain of youth" hypoth- 
esis; Perrin 2009). Because we do not yet know the 
sex-determining gene(s) of S. tropicalis or other frogs that 
might have inherited this sex determination system from a 
recent common ancestor (e.g., genera Hymenochirus, 
Pseudhymenochirus, or Pipa), we cannot determine at this 
time whether turnover or recombination best accounts for 
the apparent homomorphy of the sex chromosomes of this 
species. Additional identification of sex-determining genes in 
amphibians, and analysis of their evolutionary histories and 
genomic context, is thus an exciting direction for future 
research. 

Supplementary Material 

Supplementary figures S1-S3 and table S1 are available at 
Genome Biology and Evolution online (http://www.gbe. 
oxfordjournals.org/). 
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